{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 一.CRF的维特比算法\n",
    "\n",
    "条件随机场的预测问题是给定条件随机场$P(Y\\mid X)$和输入序列（观测序列）$x$，求条件概率最大的输出序列（标记序列）$y^*$的过程，即：   \n",
    "\n",
    "$$\n",
    "y^*=arg\\max_yP_w(y\\mid x)\\\\\n",
    "=arg\\max_y\\frac{exp(w^TF(y,x))}{Z_w(x)}\\\\\n",
    "=arg\\max_yexp(w^TF(y,x))\\\\\n",
    "=arg\\max_yw^TF(y,x)\n",
    "$$  \n",
    "\n",
    "可以发现，少了求$Z_w(x)$，计算量少了很多，上面的式子中：   \n",
    "\n",
    "$$\n",
    "w=(w_1,w_2,...,w_K)^T\\\\\n",
    "F(y,x)=(f_1(y,x),f_2(y,x),...,f_K(y,x))^T\\\\\n",
    "f_k(y,x)=\\sum_{i=1}^nf_k(y_{i-1},y_i,x,i),k=1,2,...,K\n",
    "$$  \n",
    "\n",
    "接下来将上面的$w^TF(y,x)$进行改写：   \n",
    "\n",
    "$$\n",
    "y^*=arg\\max_y\\sum_{i=1}^nw^TF_i(y_{i-1},y_i,x)\n",
    "$$   \n",
    "\n",
    "其中：  \n",
    "\n",
    "$$\n",
    "F_i(y_{i-1},y_i,x)=(f_1(y_{i-1},y_i,x,i)),f_2(y_{i-1},y_i,x,i)),...,f_K(y_{i-1},y_i,x,i)))^T\n",
    "$$  \n",
    "\n",
    "这个问题的求解，其实与之前的HMM问题一样（[链接>>>](https://nbviewer.jupyter.org/github/zhulei227/ML_Notes/blob/master/notebooks/12_10_PGM_HMM_%E9%9A%90%E9%A9%AC%E6%A8%A1%E5%9E%8B%EF%BC%9A%E9%9A%90%E7%8A%B6%E6%80%81%E9%A2%84%E6%B5%8B.ipynb)），同样是求DAG图中的最优路径问题，同样的，我们用符号$\\delta_i(l)$表示时刻$i$标签状态为$l$的所有路径中的最优值，$\\psi_i(l)$记录$i$的前一步的标签状态，下面对算法流程做一个完整的梳理\n",
    "\n",
    "#### 算法流程\n",
    "输入：模型特征向量$F(y,x)$和权值向量$w$，观测序列$x=(x_1,x_2,...,x_n)$；   \n",
    "输出：最优路径$y^*=(y_1^*,y_2^*,...,y_n^*)$  \n",
    "\n",
    ">（1）初始化   \n",
    "$$\n",
    "\\delta_1(j)=w^TF_1(y_0=start,y_1=j,x),j=1,2,...,m\n",
    "$$   \n",
    "\n",
    ">（2）对$i=2,3,...,n$   \n",
    "$$\n",
    "\\delta_i(l)=\\max_{1\\leq j\\leq m}(\\delta_{i-1}(j)+w^TF_i(y_{i-1}=j,y_i=l,x)),l=1,2,...,m\\\\\n",
    "\\psi_i(l)=arg\\max_{1\\leq j\\leq m}(\\delta_{i-1}(j)+w^TF_i(y_{i-1}=j,y_i=l,x)),l=1,2,...,m\n",
    "$$  \n",
    ">（3）终止   \n",
    "\n",
    "$$\n",
    "y_n^*=arg\\max_{1\\leq j \\leq m}\\delta_n(j)\n",
    "$$  \n",
    "\n",
    ">（4）回溯路径   \n",
    "$$\n",
    "y_i^*=\\psi_{i+1}(y_{i+1}^*),i=n-1,n-2,...,1\n",
    "$$   \n",
    "\n",
    "求得最后路径：$y^*=(y_1^*,y_2^*,...,y_n^*)$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 二.代码实现\n",
    "在上一小节的内容上追加..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.chdir('../')\n",
    "from ml_models.pgm import CRFFeatureFunction\n",
    "import numpy as np\n",
    "\n",
    "\"\"\"\n",
    "线性链CRF的实现，封装到ml_models.pgm\n",
    "\"\"\"\n",
    "\n",
    "class CRF(object):\n",
    "    def __init__(self, epochs=10, lr=1e-3, tol=1e-5, output_status_num=None, input_status_num=None, unigram_rulers=None,\n",
    "                 bigram_rulers=None):\n",
    "        \"\"\"\n",
    "        :param epochs: 迭代次数\n",
    "        :param lr: 学习率\n",
    "        :param tol:梯度更新的阈值\n",
    "        :param output_status_num:标签状态数\n",
    "        :param input_status_num:输入状态数\n",
    "        :param unigram_rulers: 状态特征规则\n",
    "        :param bigram_rulers: 状态转移规则\n",
    "        \"\"\"\n",
    "        self.epochs = epochs\n",
    "        self.lr = lr\n",
    "        self.tol = tol\n",
    "        # 为输入序列和标签状态序列添加一个头尾id\n",
    "        self.output_status_num = output_status_num + 2\n",
    "        self.input_status_num = input_status_num + 2\n",
    "        self.input_status_head_tail = [input_status_num, input_status_num + 1]\n",
    "        self.output_status_head_tail = [output_status_num, output_status_num + 1]\n",
    "        # 特征函数\n",
    "        self.FF = CRFFeatureFunction(unigram_rulers, bigram_rulers)\n",
    "        # 模型参数\n",
    "        self.w = None\n",
    "\n",
    "    def fit(self, x, y):\n",
    "        \"\"\"\n",
    "        :param x: [[...],[...],...,[...]]\n",
    "        :param y: [[...],[...],...,[...]]\n",
    "        :return\n",
    "        \"\"\"\n",
    "        # 为 x,y加头尾\n",
    "        x = [[self.input_status_head_tail[0]] + xi + [self.input_status_head_tail[1]] for xi in x]\n",
    "        y = [[self.output_status_head_tail[0]] + yi + [self.output_status_head_tail[1]] for yi in y]\n",
    "        self.FF.fit(x, y)\n",
    "        self.w = np.ones(len(self.FF.feature_funcs)) * 1e-5\n",
    "        for _ in range(0, self.epochs):\n",
    "            # 偷个懒，用随机梯度下降\n",
    "            for i in range(0, len(x)):\n",
    "                xi = x[i]\n",
    "                yi = y[i]\n",
    "                \"\"\"\n",
    "                1.求F(yi \\mid xi)以及P_w(yi \\mid xi)\n",
    "                \"\"\"\n",
    "                F_y_x = []\n",
    "                Z_x = np.ones(shape=(self.output_status_num, 1)).T\n",
    "                for j in range(1, len(xi)):\n",
    "                    F_y_x.append(self.FF.map(yi[j - 1], yi[j], xi, j))\n",
    "                    # 构建M矩阵\n",
    "                    M = np.zeros(shape=(self.output_status_num, self.output_status_num))\n",
    "                    for k in range(0, self.output_status_num):\n",
    "                        for t in range(0, self.output_status_num):\n",
    "                            M[k, t] = np.exp(np.dot(self.w, self.FF.map(k, t, xi, j)))\n",
    "                    # 前向算法求 Z(x)\n",
    "                    Z_x = Z_x.dot(M)\n",
    "                F_y_x = np.sum(F_y_x, axis=0)\n",
    "                Z_x = np.sum(Z_x)\n",
    "                # 求P_w(yi \\mid xi)\n",
    "                P_w = np.exp(np.dot(self.w, F_y_x)) / Z_x\n",
    "                \"\"\"\n",
    "                2.求梯度,并更新\n",
    "                \"\"\"\n",
    "                dw = (P_w - 1) * F_y_x\n",
    "                self.w = self.w - self.lr * dw\n",
    "                if (np.sqrt(np.dot(dw, dw) / len(dw))) < self.tol:\n",
    "                    break\n",
    "\n",
    "    def predict(self, x):\n",
    "        \"\"\"\n",
    "        维特比求解最优的y\n",
    "        :param x:[...]\n",
    "        :return:\n",
    "        \"\"\"\n",
    "        # 为x加头尾\n",
    "        x = [self.input_status_head_tail[0]] + x + [self.input_status_head_tail[1]]\n",
    "        # 初始化\n",
    "        delta = np.asarray([np.dot(self.w, self.FF.map(self.output_status_head_tail[0], j, x, 1)) for j in\n",
    "                            range(0, self.output_status_num)])\n",
    "        psi = [[0] * self.output_status_num]\n",
    "        # 递推\n",
    "        for visible_index in range(2, len(x) - 1):\n",
    "            new_delta = np.zeros_like(delta)\n",
    "            new_psi = []\n",
    "            # 当前节点\n",
    "            for i in range(0, self.output_status_num):\n",
    "                best_pre_index_i = -1\n",
    "                best_pre_index_value_i = 0\n",
    "                delta_i = 0\n",
    "                # 上一轮节点\n",
    "                for j in range(0, self.output_status_num):\n",
    "                    delta_i_j = delta[j] + np.dot(self.w, self.FF.map(j, i, x, visible_index))\n",
    "                    if delta_i_j > delta_i:\n",
    "                        delta_i = delta_i_j\n",
    "                    best_pre_index_value_i_j = delta[j] + np.dot(self.w, self.FF.map(j, i, x, visible_index))\n",
    "                    if best_pre_index_value_i_j > best_pre_index_value_i:\n",
    "                        best_pre_index_value_i = best_pre_index_value_i_j\n",
    "                        best_pre_index_i = j\n",
    "                new_delta[i] = delta_i\n",
    "                new_psi.append(best_pre_index_i)\n",
    "            delta = new_delta\n",
    "            psi.append(new_psi)\n",
    "        # 回溯\n",
    "        best_hidden_status = [np.argmax(delta)]\n",
    "        for psi_index in range(len(x) - 3, 0, -1):\n",
    "            next_status = psi[psi_index][best_hidden_status[-1]]\n",
    "            best_hidden_status.append(next_status)\n",
    "        best_hidden_status.reverse()\n",
    "        return best_hidden_status"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 2, 3, 0, 1, 3, 4]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 测试一下\n",
    "x = [\n",
    "    [1, 2, 3, 0, 1, 3, 4],\n",
    "    [1, 2, 3],\n",
    "    [0, 2, 4, 2],\n",
    "    [4, 3, 2, 1],\n",
    "    [3, 1, 1, 1, 1],\n",
    "    [2, 1, 3, 2, 1, 3, 4]\n",
    "]\n",
    "y = x\n",
    "\n",
    "crf = CRF(output_status_num=5, input_status_num=5)\n",
    "crf.fit(x, y)\n",
    "crf.predict(x[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "预测结果一样就对了..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
